/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains classes for absorption, dust models and scattering.

// $Id$


#ifndef __optical__
#define __optical__

#include "boost/thread/mutex.hpp"
#include "mcrx-types.h"
#include "ray.h"
#include "random.h"
#include "constants.h"
#include "angular.h"
#include "mono_poly_abstract.h"
#include "mcrx-debug.h"

namespace mcrx {
  template <typename, typename> class optical_data;
  template <typename> class absorber;
  template <typename> class intensity;
  template <typename, typename> class intensity_pointer;
}

/** Container class for the grid cells. This class contains an emitter
    and an absorber as template classes, in case the grid_cells need
    to contain both absorbers and emitters. */
template <typename emitter_type, typename absorber_type>
class mcrx::optical_data {
public:
  typedef emitter_type T_emitter;
  typedef absorber_type T_absorber;  
private:
  T_emitter em;
  T_absorber ab;

public:
  optical_data (const T_emitter& e, const T_absorber& a):
    em (e), ab (a) {};

  // these do not return const references, I think we might need to
  // assign to them...
  T_emitter& get_emitter () {return em;};
  T_absorber& get_absorber () {return ab;};

  //  copy constructor and assignment operator are implicitly generated

  // This is needed for the grid_factory
  optical_data& operator+= (const optical_data& rhs) {
    em+= rhs.em; ab+= rhs.ab; return *this;};
  optical_data& operator*= (const optical_data& rhs) {
    em*= rhs.em; ab*= rhs.ab; return *this;};
  optical_data operator*(const optical_data& rhs) {
    return optical_data (*this)*= rhs;}; 
  static optical_data unification (const optical_data& sum, int n) {
    return optical_data (T_emitter::unification (sum.em, n),
			 T_absorber::unification (sum.ab, n));};
};


/** Describes the absorbing/scattering medium in a grid cell. It
    contains an array of densities of the different scattering
    species, and the intensity of the radiation field in the cell. */
template <typename lambda_type>
class mcrx::absorber {
public:
  typedef mcrx::T_float T_float;
  typedef lambda_type T_lambda;
  typedef ray<T_lambda> T_ray; 
private:
  /** Density vector contains the densities of the scatterers that are
      specified in a dust_model object. The grid itself has no
      knowledge of the dust_model used, it is kept by the xfer
      object. */
  T_densities rho_;

  /** Momentum density of medium (rho*v), for kinematic
      calculations. No velocity segregation between the different
      species is possible. */
  vec3d rhov_;

#ifdef MCRX_DEBUG_LEVEL
  /// ray contents absorbed. Only used when debugging.
  T_lambda absorbed_; 
#endif

  mutable boost::mutex absorption_mutex;

  /** Mean radiation intensity, as calculated from the rays traversing
      the cell. (Actually, it's sum of the extensive quantity dl*I for
      all rays that have traversed the cell. The intensity is
      calculated from this by sum(I*dl)/4pi*V.) */
  T_lambda intensity_;

public:
  /** Constructor takes the density array and also an optional zero
      value for the intensity and absorption members.  This only makes
      sense if they are arrays, because they need to be sized.  We are
      not using independent_copy for the intensity, because we might
      be feeding weak-referenced arrays that should be copied by
      reference.  */
  absorber (const T_densities& rho, const vec3d& v, 
	    const T_lambda& intensityzero=T_lambda()):
    rho_(rho), rhov_(v*sum(rho)),
#ifdef MCRX_DEBUG_LEVEL
    // absorbed needs to be independent from intensity
    absorbed_(independent_copy(intensityzero)),
#endif
    intensity_(intensityzero)
  {
    // bogus?
    // Note that we DO NOT want to use "independent_copy" here because
    // it would make it impossible to create absorbers referring to a
    // master array of intensities and absorption.
    assert (condition_all(rho_ >= 0));

    // however, we do want to set the intensity to zero, though...
    intensity_ = 0;
#ifdef MCRX_DEBUG_LEVEL
    absorbed_ = 0;
#endif
  };
 
  /** We need to write the copy constructor explicitly because we need
      to be sure we don't make reerence copies if T_lambda is an arry
      \todo IS THIS TRUE? What if they are pointing to a master array?
      Maybe disallowing copying would be better?  */
  absorber (const absorber& rhs): 
    rho_(independent_copy(rhs.rho_)), rhov_(rhs.rhov_),
#ifdef MCRX_DEBUG_LEVEL
    // absorbed needs to be independent from intensity
    absorbed_(independent_copy(rhs.absorbed_)),
#endif
    intensity_(independent_copy(rhs.intensity_)) {};

  /** This accessor is so that a grid can use absorbers both directly
      and with the the optical_data class (which is now unnecessary in
      most cases), while still retaining that possibility. Hopefully
      use of absorber gets optimized out...  */
  absorber& get_absorber() {return *this;};
  const absorber& get_absorber() const {return *this;};

  /// Return densities.
  const T_densities& densities () const {return rho_;};
  /// Return velocity.
  const vec3d velocity () const {return rhov_/sum(rho_);};
#ifdef MCRX_DEBUG_LEVEL
  /// Return absorbed luminosity, if debugging.
  const T_lambda& absorbed () const {return absorbed_;};
#endif
  /// Return intensity vector.
  const T_lambda& intensity () const {return intensity_;};

  /** Add to the intensity array. This function is protected by the
      mutex. */
  template <typename T>
  /** Add to the absorbed array. This function is protected by the
      mutex. */
  void add_intensity (const blitz::ETBase<T>&);
  void add_absorbed (const T_lambda&);
  /** Add to the absorbed array. This function is protected by the
      mutex. */
  template <typename T>
  void add_absorbed (blitz::_bz_ArrayExpr<T>);

  /** Sets the absorption member.  If T_lambda is an array, it also
      resizes it appropriately.  Note that this observation is not
      protected by the mutex since it is rarely used.  */
  void set_absorbed (const T_lambda& a) {
    DEBUG(1,assign (absorbed_, a););};

  /** Sets the intensity member.  If T_lambda is an array, it also
      resizes it appropriately.  Note that this observation is not
      protected by the mutex, since it is rarely used.  */
  void set_intensity (const T_lambda& i) {assign (intensity_, i);};
  
  /// \name Methods necessary for the grid_factory.
  ///@{
  absorber& operator+= (const absorber& rhs);
  absorber& operator*= (const absorber& rhs) {
    // element-wise multiplication of all quantities
    rho_ *= rhs.rho_;
    rhov_ *= rhs.rhov_;
    DEBUG(1,absorbed_*= rhs.absorbed_;);
    intensity_*= rhs.intensity_;
    return *this;};
  absorber operator*(const absorber& rhs) {
    return absorber (*this)*= rhs;};
  static absorber unification (const absorber& sum, int n) {
    absorber a (sum);
    // rho is intensive
    a.rho_ /= n;
    // rho*v is also intensive
    a.rhov_ /= n;

    // all others are extensive, so the sum is correct
    return a;};
  ///@}
};


// ***function definitions ***

template <typename ray_content_type>
mcrx::absorber<ray_content_type>&
mcrx::absorber<ray_content_type> ::operator+= (const absorber& rhs)
{
  rho_ += rhs.rho_;
  rhov_ += rhs.rhov_;
  intensity_ += rhs.intensity_;
  DEBUG(1,absorbed_ += rhs.absorbed_;);
  return *this;
}


/** Abstracts the data type used to represent intensity in
  the grid cells. not used we just use straight T_lambda */
template <typename intensity_type>
class mcrx::intensity {
  intensity_type i;
public:
  intensity (): i (0) {};
  intensity (intensity_type ii): i (ii) {};
  intensity_type value () const {return i;};
  intensity& operator+= (intensity rhs) {
    i+= rhs.i; return *this;}
};


/// Add to the radiation intensity in the absorber.
template <typename ray_content_type>
template <typename T>
void 
mcrx::absorber<ray_content_type>::
add_intensity (const blitz::ETBase<T>& i) 
{
  const T& intensity (i.unwrap());

  boost::mutex::scoped_lock absorption_lock (absorption_mutex);
  if(!same_size (intensity_, intensity)) {
    resize_like (intensity_, intensity);
    intensity_=0;
  }
  intensity_ += intensity;
}

template <typename ray_content_type>
void 
mcrx::absorber<ray_content_type>::
add_absorbed (const T_lambda& a) 
{
#ifdef MCRX_DEBUG_LEVEL
  boost::mutex::scoped_lock absorption_lock (absorption_mutex);
  if(!absorbed_.size()) {
    absorbed_.resize(a.shape());
    absorbed_=0;
  }
  absorbed_ += a;
  // intensity may be negative here for perturbative ir runs
#endif
}


template <typename ray_content_type>
template <typename T>
void 
mcrx::absorber<ray_content_type>::
add_absorbed (blitz::_bz_ArrayExpr<T> a) 
{
#ifdef MCRX_DEBUG_LEVEL
  boost::mutex::scoped_lock absorption_lock (absorption_mutex);
  if(absorbed_.size())
    absorbed_ += a;
  else {
    const T_lambda temp(a);
    absorbed_.resize(temp.shape());
    absorbed_ += temp;
  }
  assert(condition_all(absorbed_>=0));
#endif
}
  

#endif


